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Abstract 

A new fully implicit, time accurate algorithm suitable 
for chemically reacting, viscous flows in the transonic-to- 
hypersonic regime is described. The method is based on 
a class of Total Variation Diminishing (TVD) schemes 
and uses successive Gauss-Seidel relaxation sweeps. The 
inversion of large matrices is avoided by partitioning 
the system into reacting and nonreacting parts, but still 
maintaining a fully coupled interaction. As a result, the 
matrices that have to be inverted are of the same size 
as those obtained with the commonly used point implicit 
methods. In this paper we illustrate the applicability of 
the new algorithm to hypervelocity unsteady combustion 
applications. We present a series of numerical simula- 
tions of the periodic combustion instabilities observed 
in ballistic-range experiments of blunt projectiles flying 
at subdejtonative speeds through hydrogen-air mixtures. 
The computed frequencies of oscillation are in excellent 
agreement with experimental data. 

Introduction 

The development of hypersonic propulsion systems 
requires efficient computational fluid dynamics (CFD) 
codes to support experimental tests and extrapolate the 
results to full scale and flight conditions (where static 
pressures may be an order of magnitude higher than those 
achievable in experimental facilities). The acquisition of 
data in hypervelocity facilities is difficult because of the 
short test times available. In a pulse facility, for example, 
the duration of the transient test flows is on the order 


of one millisecond, which may be insufficient for the es- 
tablishment of the viscous dominated fluid dynamic and 
chemical reaction processes 1 . In a ram accelerator 2 , there 
is considerable interest in studying the transient develop- 
ment of the shock-induced combustion process, its stabil- 
ity and interactions with a reacting or nonreacting bound- 
ary layer. Due to the highly stretched grids necessary to 
resolve boundary layers and the high flow velocities in- 
volved, simulation of these problems requires a fully im- 
plicit time accurate CFD code in order to maintain high 
computational efficiency. Recent time accurate simula- 
tions of hypersonic reacting flows 3-5 have been conducted 
predominantly with explicit or point implicit methods, in 
which only the source term is treated implicitly. The rea- 
son is partly due to the fact that the Jacobian matrix 
associated with the chemistry source terms complicates 
the solution procedures for the commonly used implicit 
ADI methods because the governing equations become 
stiff. Since explicit or point implicit methods are con- 
strained by the CFL condition, they are very inefficient 
for solving viscous reacting flows. 

In the present work, we have developed a new relax- 
ation method that is used with an implicit TVD scheme 
to solve the fully coupled chemically reacting flow equa- 
tions. The algorithm is based on successive backward and 
forward Gauss-Seidel relaxation sweeps. The inversion of 
large matrices is avoided by partitioning the system into 
reacting and nonreacting parts, but a fully coupled in- 
teraction is still maintained. As a result, the matrices 
that have to be inverted are of the same size as those 
obtained with point implicit methods. The method pre- 
sented is second order accurate in time and space, and 
allows variable time stepping. The time step is automat- 
ically selected based on the satisfaction of a convergence 
criterion. 

The performance of the new method is illustrated by 
applying it to hypervelocity, unsteady, shock-induced 
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combusting flows. We present a series of numerical sim- 
ulations of the periodic combustion instabilities observed 
in ballistic-range experiments of blunt projectiles flying 
at subdetonative speeds through hydrogen-air mixtures. 

Numerical Formulation 
Governing Equations 


then a second-order, variable-step backward differentia- 
tion formula (BDF) applied to Eq. 1 gives 

AQJ, = 7 AQ£ 1 -/JAfTF, + , Jt 

( 4 ) 

where the BDF method coefficients 7 and /? are 


We consider the nonequilibrium Navier-Stokes equa- 
tions for two-dimensional or axisymmetric flow, in which 
the global continuity equation is replaced by species con- 
tinuity equations. They can be expressed in the following 
conservation form for a gas mixture containing n species 
and in general curvilinear coordinates (£,1?) 

s -s.)=w a. 



A 2 

1 + 2A 
(f + A) 
1 + 2A 


( 5 ) 

( 6 ) 


The function F J+ | fc is the numerical flux in the £ di- 
rection evaluated at (j + 1 , k). It is computed using Yee’s 
second order TVD scheme 7 as follows. Let A = <9F/dQ 
and B = dG/dQ. Then, for a pseudo finite volume for- 
mulation, the numerical flux can be expressed as 


where 


[SI 


Q = J- 1 
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(2) 


is the vector of dependent variables. 

The equations describe two-dimensional flow if j = 0 
and axisymmetric flow if j = 1. The variables are the 
mass density of the tth species Pi (mixture mass density 
p = P*)i velocity components u and v and the 
total energy per unit volume e. J denotes the grid Ja- 
cobian and F and G are the inviscid flux vectors in the 
£ and 17 directions respectively. Similarly, F„ and G„ 
are the viscous fluxes. The terms S and S„ are the ax- 
isymmetric source terms, and W is the chemical source 
term. A detailed description of all the terms appearing in 
Eq. ( 1 ), and of the additional state and constitutive equa- 
tions needed to close the system are given by Yungster 6 . 


Description of the Implicit TVD Scheme 


+ F i+i,fc + R i+i*j+i) ( 7 ) 

In two dimensions, this pseudo finite volume formulation 
can be made “truly” finite volume by a slight modification 
in the calculation of some of the metric terms. See Ref. 7 
for further details. 

Here, R J+ j denotes the matrix of eigenvectors of the 
flux Jacobian matrix A evaluated at some symmetric av- 
erage of Qj,fc and Qi+i,*, denoted as Q J+ ^. Abgrall’s ex- 
tension of Roe’s averaging for real gas mixtures 8 is used 
in the present study. Similarly, one can define in this 
manner the numerical flux in the i] direction. 

The elements, of the dissipation vector are: 

05 +} = -^ + j)K + j-O‘- +i ] (8) 

Here a'. + j denotes the eigenvalues of A evaluated at 
Qj + j , and a‘. + ^ denotes the elements of the vector <*j+y- 

®' ; +j(Qj+l.fc* 7 i+l,fc — Q},kJj,k) 

a,+i = 0.5(Jj+i,fc + Jj^k) (9) 

The function ’F is : 


For simplicity and without loss of generality, the im- 
plicit method is presented here for the two-dimensional 
Euler equations (F„ = G„ = 0 ). For the Navier-Stokes 
equations, a central difference approximation for the vis- 
cous fluxes should be added to the right hand side op- 
erator, and the corresponding viscous Jacobians should 
be added to the implicit operator. Let the time step 
At" = t n+1 - t", and the change AQ" = Q" +1 - Q", 


and define 


A = 


At" 

At"- 1 


( 3 ) 
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The entropy correction parameter, e, in Eq. ( 10 ) is taken 
to be a function of the velocity and sound speed 7 


e j+] — * 


4[(vlH¥+v^i +) J u > 
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where a is the frozen sound speed, U and V are the con- 
travariant velocities (U — £*u+£ y u, V = and e 

is a small number in the range 0 < e < 0.4, which controls 
the convergence rate and the sharpness of discontinuities. 

The “limiter” function Qj+^ used in this study is the 
following: 

Q l j+ i = 2 a l j+i , 2a} +} , + a* + j)] 

( 12 ) 

The eigenvalues and eigenvectors of the fully coupled 
chemically reacting equations, as well as the resulting ex- 
pressions for R# are given in Ref. 9. 

Iterative Scheme 


We can rearrange Eq. 13 by defining “plus” and “minus” 
matrices having nonnegative and nonpositive eigenvalues 
respectively, 

A+ - 2* A,-1,t + 

(23) 


(24) 

A ~ = 2 ■ A J +1 * - 

(25) 

B- = ~(B M .,-n’ i+} ) 

(26) 

Define the matrix D as 



Following the work of Yee 10 , we linearize Eq. 4 in a 
conservative manner: 


- C,, t ]"E> = -RHS (13) 

where RHS is the right-hand side of Eq. 4, with the ex- 
ception that the term in square brackets is evaluated at 
time level n instead of n + 1. Also 

E n = Q n+1 - Q n (14) 


and 



(15) 


The remaining terms in Eq. 13 are defined as follows: 



[A j+1 , k - 


( 16 ) 


f3At n D = I+/3At«(Of +iifc + ■ +Oj^ 4 ■ + «?*_*) 

(27) 

Then Eq. 13 reduces to the following expression: 

- B + E" t _ 1 + (D — C)E" it + 

A_ E" + i,t + B-E” >+1 = -^ (28) 

This equation is solved iteratively using a lower-upper 
relaxation procedure consisting of successive Gauss-Seidel 
(LU-SGS) sweeps as follows: 

Step 1 


(D-C)E;,, = A-E'.^+B+EJj,.,-^- 


followed by Step 2 



(17) 

where ft is defined as 


= ^+l A y+i^,'+i 

(18) 


(19) 

(Here R fc+ ^ denotes the matrix of eigenvectors of the flux 
Jacobian matrix B evaluated at some symmetric average 
of Qj,k and Qj.t+i). The terms A^ + , and A^ are 
diagonal matrices defined as 

A^ +i =<fio<7[^(o( +i )] 

(20) 

Afc + i = 

(21) 


The nonstandard notation of Yee 7 is used. That is 


Hf +iifc E" = i[A i+1)fc E? +M - d] +hk (E] +hk - E" jfc )] 

( 22 ) 


(D-C)(E;.J-E* fc ) = -[A-(E<^ t -E^>) + 

B-(e5“Vi-EK ) )1 (30) 

Here, (m) is the subiteration counter, and we have used 
the notation 

E (m) s (E n ) (m) (31) 

Note that Eqs. 29 and 30 require only block diago- 
nal matrix inversions — a significant advantage over ADI 
based methods, which require the solution of block tridi- 
agonal or block pentadiagonal systems. To avoid the 
inversion of large matrices, we can modify the iterative 
scheme given in Eqs. 29 and 30 by partitioning the vector 
Q and the matrices D and C as follows: 


Q 7 1 


■ z 11 z 12 ‘ 

Q" J 
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Here, Q 7 is a vector of length n r , where 7V is the num- 
ber of reacting species that participate in the chemical 



reactions. The vector Q 11 includes all the remaining vari- 
ables. Then, the LU-SGS method given by Eqs. 29 and 
30 can be simplified as 
Step 1 


z 22 (E* fc )" = [a+e^+b+e;,*.,- 


RHS 

0At n 


1) 


A"E 




B — 

I “ 


5 2 1 ( E (m- 1) ) / 


(33) 


RHS 


Z U (E U) 1 = [A + E*_ 1)fc + B+E^.j - ^ Afn 

a-eJ^-b-eJ^V- 

Z 12 (E£ fc )" (34) 


followed by Step 2 

z 22 (e<J-e;,*)" = -[A-(E5> ir Eg>) + 

B"(E^i-e5^)]"- 

Z^E^-Eg- 1 *)' (35) 

Z 11 (E(J-E;, fc ) f = -[A-(E^ tJfe -E^ ) ) + 
b^e^-e^)]'- 

Z»(Eg ) -EJ tfc )" (36) 

Therefore, at every grid point and for steps 1 and 2, we 
solve first a (»+3-n P ) x (n+3-»r) system of equations 
for the “nonreacting” variables (i.e., Eq. 33 in step 1), 
followed by a solution of a (nr) x (nr) system for the 
reacting species (Eq. 34), and using the most updated 
values available for the nonreacting variables. Because 
the operation count is proportional to N 3 when solving 
an N x N system of equations, this modified iterative 
scheme is preferred over the unpartitioned scheme given 
by Eqs. 29 and 30, which requires the solution of (n + 
3) x (n + 3) systems of equations. 

The solution of the block system of equations is done by 
LU-decomposition of the Z n and Z 22 matrices followed 
by forward substitution and back substitution. The LU- 
decomposition is stored at every grid point so that at 
subsequent iterations only the forward and back substi- 
tutions are required. 

The iterative process given by Eqs. 33-36 is repeated 
until convergence is achieved. The procedure can be vec- 
torized by taking sweeps along diagonal directions. That 
is, for a £ - 77 rectangular computational domain, step 1 
above sweeps the domain from the lower left corner to 
the upper right comer, and step 2 takes sweeps from the 
upper right comer back to the lower left comer of the 
computational domain. 


Convergence Test 


The test for iteration convergence is similar to that used 
in the LSODE code 11 , a solver for initial value problems in 
ordinary differential equation systems. The test is based 
on the successive differences (Q< ro ) - Q (m-1 >), as com- 
pared with an error weight vector, EWTi, defined below. 
Convergence is said to occur if 



at every grid point, where 


EWT i = RTOL i \Q < i m) \ + ATOLi-, t = (38) 


and BTOLi and ATOLi are, respectively, the user- 
supplied local relative and local absolute error tolerances 
for the i th component, and Cl is a factor of order one. If 
convergence is achieved in less than the maximum num- 
ber of iterations permitted, (iter ) ma *, then the time step 
is increase by 5% — 10% for the next integration step. If 
the iterative scheme fails to converge in (iter) mox , then 
the time step is reduced by 20% — 50%. (Our experience 
to date suggests use of (iter) max — 5.) 


Numerical Simulations 


Description of test problem 

The performance of the new algorithm is demonstrated 
by a series of numerical simulations of spherical-nosed 
projectiles fired into hydrogen-air mixtures at subdeto- 
native speeds. Experimental studies of such flows were 
conducted in the 1960’s and early 1970’s, most notably 
by Ruegg and Dorsey 12 , Behrens et al. 13 and Lehr 14 ’ 15 . 
These experiments indicated that the flow field is highly 
unstable under a wide range of conditions. Depending 
on the flow conditions, the periodic instability appears 
in two main types referred to as “regular regime” and 
“large-disturbance” regime. In the present paper, we fo- 
cus on the regular regime characterized by a highly regu- 
lar, periodic flow structure (see Figs. 3a, 5a and 7a, taken 
from the work of Lehr 14 ). A mechanism for the regular 
regime periodic combustion instability, containing most 
of the essential features, was first proposed by Toong and 
coworkers 16 ’ 17 . Three important processes form the ba- 
sis of Toong’s model. They are described below, and 
are schematically illustrated in the one-dimensional di- 
agrams shown in Fig. 1. This one-dimensional analysis 
could be viewed as representing the flow along the stag- 
nation streamline of the blunt projectile. 

(I) A new reaction front generates compression waves, 
which have been termed “reaction shocks”, that travel 
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upstream (towards the bow shock) and downstream (to- 
wards the body stagnation point) as shown in Fig. la. 
Alpert and Toong 17 demonstrated that such waves are 
necessary for mass and momentum conservation. The 
upstream facing reaction shock eventually must overtake 
the bow shock, strengthening it and producing a reflected 
rarefaction wave and a contact discontinuity. The rarefac- 
tion wave is very weak and plays no role in the instability 
mechanism. It is usually neglected. 

(II) Across the contact discontinuity, which marks the 
boundary of the gas which has passed through the rein- 
forced bow shock, there is a sharp increase in tempera- 
ture. The induction time in the hot gas is shortened and 
as a result, a new reaction front is created between the 
bow shock and the original reaction front, as shown in 
Fig lb. This new reaction front generates a pair of reac- 
tion shocks as discussed in (I). The reacted gas behind the 
new reaction front eventually reaches the location of the 
original reaction front. The latter is extinguished because 
no more unreacted gas is present. When a reaction front 
is extinguished, it must be accompanied by the genera- 
tion of upstream and downstream rarefaction waves which 
have a strength comparable to the reaction shocks 17 . 

(III) The rarefaction wave generated by the extinguishing 
of the original reaction front eventually overtakes the bow 
shock. As the rarefaction wave penetrates the bow shock, 
a fan of weak pressure waves are reflected and a zone of 
decreasing temperature is created (Fig. lc). As a result, 
the ignition delay progressively increases and the reaction 
front recedes towards the body. 

The above processes form the basis of the instability 
mechanism. The complete cycle can be explained, based 
on these processes, in a x-t diagram as shown in Fig- 
ure 2. This figure shows the McVey-Toong mechanism 16 , 
with the addition of the reflected shocks from the body. 
The importance of these reflected shocks has been recog- 
nized only recently, following the computational studies 
of Wilson and Sussman 3 , and Matsuo et a/. 4 . The x-t 
diagram of Fig. 2 shows a particular case for which the 
reflected shock from the body overtakes the bow shock 
at exactly the same time as the reaction shock does (res- 
onant mode). In general, this will not be the case and 
the diagram will become more complex. The cycle of 
events can be described as follows: A reaction shock and 
a reflected shock overtake the bow shock simultaneously, 
strengthening it and reflecting a contact discontinuity and 
a weak rarefaction wave not shown. A zone of decreasing 
temperature follows the contact discontinuity, due to a 
rarefaction wave penetrating the bow shock. The origin 
of this rarefaction wave is explained below. The hot gas 
behind the contact discontinuity begins to react creating 
a new reaction front and a pair of reaction shocks moving 
upstream (towards the bow shock), and downstream (to- 
wards the body). When the contact discontinuity reaches 


the original reaction front at a later time, it extinguishes it 
and generates a pair of rarefaction waves. The new reac- 
tion front gradually recedes due to the zone of decreasing 
temperature until it reaches its original location. When 
the reaction shock (and the reflected shock) overtake the 
bow shock, a new cycle will begin. 

The computation of the periodic combustion instabili- 
ties in exothermic blunt body flows represents an excel- 
lent test problem, because the numerical method must be 
able to accurately resolve all the complex wave interac- 
tions previously described. The numerical simulation of 
ballistic range experiments of projectiles fired into det- 
onable mixtures can be traced to the work of Yungster 
et a/. 18 ’ 19 , Lee & Deiwert 20 , Wilson & MacCormack 21 , 
and Soetrisno & Imlay 22 . These studies focused on the 
cases that exhibited stable combustion. Numerical sim- 
ulations of the pulsating shock-induced combustion flows 
were first reported by Wilson and Sussman 3 . Subsequent 
studies have been conducted by Matsuo et at.*, Ahuja & 
Tiwari 5 and by Hosangadi et al. 23 . There are substan- 
tial differences in the results predicted by these authors 
(such as frequency and amplitude of the oscillations) even 
though all studies were based on similar combustion mod- 
els. These differences indicate that the problem may be 
quite sensitive not only to the reaction mechanism used, 
but also to the numerics and grid used. The work of Wil- 
son and Sussman was based on a flux-vector splitting or 
upwind TVD point implicit scheme (i.e., the convective 
terms were treated explicitly whereas the source term was 
treated implicitly), and utilized the logarithmic transfor- 
mation of the species conservation equations. Matsuo et 
al. utilized a point implicit TVD scheme, Ahuja & Ti- 
wari utilized the SPARK2D code based on MacCormack’s 
scheme with added artificial dissipation. The study of 
Hosangadi et al. was conducted using the CRAFT code 
which has an implicit Roe/TVD based upwind formula- 
tion with ADI factorization. 

Results and Discussion 

A numerical simulation of Lehr’s 14,15 ballistic range 
experiments is presented. These consisted of spherical- 
nosed projectiles having a diameter of 15-mm fired into 
a premixed, stoichiometric mixture of hydrogen-air. The 
pressure and temperature of the gas mixture was 320 torr 
and 293 K respectively for all the cases considered. The 
Chapman-Jouguet detonation Mach number for this mix- 
ture is Mcj = 5.10. The combustion model used in the 
present study is based on the Jachimowski 24 mechanism 
excluding the nitrogen reactions. This mechanism is iden- 
tical to the one used by Matsuo et al.*, and is given in 
Table 1. The calculations were carried out on a 220 x 220 
uniform grid, and viscous effects were neglected. The 
periodic instabilities have been experimentally observed 
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only for subdetonative speeds, i.e., for Mach numbers be- 
low Mcj. Computational results are presented for three 
Mach number flows, M = 4.18,4.48 and 4.79. 

Figure 3 shows instantaneous density, pressure and 
water mass fraction contours for the M = 4.18 case. 
For comparison, the experimental shadowgraph image of 
Lehr 15 is also shown in Fig. 3. The water mass fraction 
contour plot (Fig. 3d) indicates the boundary of the reac- 
tion front. The density contour shows the relatively large 
induction zone separating the bow shock and the combus- 
tion front. The large induction time in this case, produces 
a low pulsating frequency with only one visible pulse in 
the nose region of the projectile. This prediction agrees 
with the experimental result shown in Fig. 3a. Figure 3b 
also shows clearly the contact discontinuity emanating 
from the edge of the reaction front pulse and moving to- 
wards the bow shock. The fact that this wave constitutes 
a contact discontinuity can be deduced by its absence 
from the pressure contour plot (Fig. 3c). 

The history of density and pressure along the stagna- 
tion streamline can be plotted in an x - 1 diagram corre- 
sponding to that in Fig 2. The computed x - t diagram 
for the M = 4.18 case is presented in Fig 4. This fig- 
ure shows two pulses in the reaction front. The creation 
of new reaction fronts as well as reaction shocks and the 
contact discontinuity can be clearly identified in the den- 
sity plot. This plot shows also the existence of multiple, 
weak waves reflecting from the projectile surface and bow 
shock. The presense of these weak waves can be explained 
by the lack of synchronization between the reaction shock 
moving towards the bow shock and the reflected shock 
from the projectile surface. 

Table 2 shows the frequency of oscillation reported in 
the experimental work of Lehr 15 , compared with the re- 
sults obtained in the present study and previous com- 
putations. The frequency of oscillation predicted in the 
present work for the M = 4.18 case is 163 KHz, which is 
in excellent agreement with both the experimental value 
of 148 KHz 15 and the computational result of Matsuo et 
al. 4 , who predicted a frequency of 160 KHz. 

At higher Mach number flows, the frequency of oscilla- 
tion will increase, because, as seen from Fig 2, the period 
of one oscillation is the time required for a reaction shock 
to overtake the bow shock plus the time required for the 
contact discontinuity to convert from the bow shock to 
the new reaction front (i.e. the induction time). The 
higher the Mach number is, the higher the temperature 
behind the bow shock, and the shorter the induction time. 

Figure 5 shows instantaneous density, pressure, and wa- 
ter mass fraction contours, and the corresponding shad- 
owgraph experimental image for a Mach number flow of 
M = 4.48. The shadowgraph picture for this case was 
obtained with a camera located at a 12° angle from the 
perpendicular direction. Four pulses can be clearly iden- 


tified. A contact discontinuity emanating from each of 
the pulses is clearly visible in the density plot (Fig. 5b). 

The pressure plot (Fig. 5c) shows also a pair of reaction 
shocks, one moving towards the bow shock and the other 
towards the projectile surface. The fact that no reaction 
shocks are observed downstream indicates that the insta- 
bility mechanism is driven primarily by the interactions 
occurring in a small region near the stagnation point. The * 
H 2 O mass fraction contour plot in Fig. 5d confirms the 
formation of the new reaction fronts. 

A typical calculation required approximately 4000 steps * 
and 7.5 hrs of CPU time on a CRAY-C90. The number 
of subiterations ( iter) max was 5, and the maximum CFL 
number varied from 2 to 5 approximately. 

The history of density and pressure along the stagna- 
tion streamline is plotted in Fig. 6. Note that for this 
condition the wave interaction is almost perfectly syn- 
chronized (resonant mode), that is, the reflected shock 
from the body overtakes the bow shock almost at the 
same time as the reaction shock. This situation corre- 
sponds exactly to the schematic wave diagram of Fig. 2. 

The near perfect synchronization between the reflected 
shock and the reaction shock produces a very clean wave 
pattern of nearly constant frequency and amplitude. The 
frequency predicted by the present calculation is 431 KHz 
and the flow is exactly periodic. This is in excellent agree- 
ment to the reported experimental frequency of oscillation 
of 425 KHz 15 . 

The final calculation considered a Mach number of 
M = 4.79. Figure 7 presents instantaneous density, pres- 
sure, and water mass fraction contours, as well as the ex- 
perimental shadowgraph image (obtained from a perpen- 
dicular direction as in Fig. 3a). The reaction front is much 
closer to the bow shock, and therefore the frequency of the 
oscillation is higher. Figure 8 shows the history of density 
and pressure along the stagnation streamline. Note that 
in this case the reflected shock overtakes the bow shock 
after two periods of oscillation (second resonant mode). 
However, the reaction shock and the reflected shock are 
not perfectly synchronized as in the previous case. As 
a result, the wave structure is not as clean as for the 
M = 4.48 case, and the frequency and amplitude of the 
pulses is not exactly constant. The predicted frequency 
of oscillation varied from 701 to 716 KHz. This result is 
also in very good agreement with the experimental fre- 
quency of 720 KHz 15 , and with the computational result 
of Matsuo et al. 4 , who predicted a frequency of 725 KHz. 

The computation of Matsuo et al. predicted an exactly [ 
second resonant mode. This case was also computed by 
Wilson and Sussman 3 and Hosangadi et al. 23 . The former 
predicted a frequency of oscillation of 530 KHz, and the • 
latter a frequency in the 450-500 KHz range, both sig- 
nificantly lower than those obtained in the present study 
and by Matsuo et al. and Lehr 15 . One possible cause 
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for the discrepancy in the results of Hosangadi et al., is 
the reaction mechanism. Their study utilized a simpler 
7-spedes model that may not be accurate enough. The 
study of Wilson and Sussman included 13 nitrogen reac- 
tions in addition to the reactions considered in the present 
study (and in Matsuo’s work). To test the effects of the 
nitrogen reactions, we conducted static ignition studies 
of stoichiometric hydrogen-air mixtures at the postshock 
conditions, with and without the nitrogen reactions con- 
sidered by Wilson and Sussman. Our results showed that, 
for the conditions considered, the nitrogen reactions have 
no effect on the ignition delay, rate of heat release or 
equilibrium temperature. Therefore, the difference in the 
results cannot be attributed to the reaction model. 

Finally, it is interesting to compare the stagnation pres- 
sure variation as a function of time for the three cases 
studied. This comparison is shown in Fig. 9. Note that 
for the nearly resonant M = 4.48 case the pressure pat- 
tern repeats itself almost identically. Also, the amplitude 
of the oscillation is much larger than the other two cases 
as a result of the reinforcing, synchronized action of the 
reaction and reflected shocks. In the M = 4.79 case one 
can notice the change in the pressure pattern as a result of 
the lack of synchronization. Also, there is a larger vari- 
ation in the amplitude of the oscillations. The average 
pressure level increases with Mach number as expected. 

Conclusions 

A new time-accurate algorithm suitable for chemically 
reacting flows in the transonic-to-hypersonic regime was 
presented. The salient features of the new method are: 
1) fully implicit, fully coupled scheme; 2) second-order 
accurate in time and space; 3) automatic time-step selec- 
tion; and 4) an efficient iterative scheme that avoids the 
inversion of large matrices. 

Numerical simulations of periodic combustion instabil- 
ities observed in ballistic range experiments indicate that 
the new method can accurately and efficiently solve com- 
plex, unsteady, chemically reacting flows. 
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Table 1: U 2 -air reaction mechanism 0 



No. 

Reaction 

A 


e 

6 . 

1 

H 2 4- O 2 ^ H0 2 + h 

1.00 X 

To 1 * 

28197.38 

0.0 

2 

H 4~ O 2 ^ OH -f O 

2.60 x 

10“ 

8459.21 

0.0 

3 

H 2 + 0 ^ OH + H 

1.80 x 

10 10 

4481.37 

lsft 

4 

H 2 4- OH ^ H -+■ S 2 O 

2.20 X 

10 13 

2593.15 

0.0 

5 

OH 4* OH ^ (!) + H 2 O 

6.30 x 

10 12 

548.84 

0.0 

6 

H + OH + M^ H 2 0 + M 

2.20 X 

10 22 

0.0 

-2.0 

7 

h + h + m*±h 2 + m 

6.40 x 

10 17 

0.0 

-1.0 

8 

H + 0 + M*±0H + M 

6.00 x 

10 16 

0.0 

-0.6 

9 

H -f - O 2 + Af ^ HO 2 4* Af 

2.10 x 

10 15 

-503.52 

0.0 

10 

O + 0 + M O 2 + M 

6.00 x 

10 13 

-906.34 

0.0 

11 

HO 2 + H?±OH + OH 

1.40 x 

10 14 

543.81 

0.0 

12 

H0 2 4- H ^ H 2 0 4- O 

1.00 x 

10 13 

543.81 

0.0 

13 

HO 2 4- O ^ 0 2 4- OH 

1.50 x 

10 13 

478.35 

0.0 

14 

H 0 2 + OH ^ H 2 0 4 ' 0 2 

8.00 x 

10 12 

0.0 

0.0 

15 

HO2 + HO2 ^ H2O2 + o 2 

2.00 x 

10 12 

0.0 

0.0 

16 

H 4- H2O2 ^ H2 4* H 0 2 

1.40 x 

10 12 

1812.69 

0.0 

17 

O 4- H2O2 ^ Off 4- HO2 

1.40 x 

10 13 

3222.56 

0.0 

18 

OH 4~ H2O2 ^ H 2 0 + HO2 

6.10 x 

10 12 

720.04 

0.0 

19 

H 2 O 2 4~ M ^ OH 4“ OH 4 - M 

1.20 x 

10 17 

22910.37 

0.0 


“Forward rate coefficient kf — i47 ( exp(— 6/7); units are moles, 
seconds, centimeters, and Kelvins. 

Third-body efficiencies: 

(6) fT 2 0 = 6.0 

(7) B 2 0 = 6.0, H 2 = 2.0 

(8) H t O = 5.0 

(9) BtO = 16.0, Hi = 2.0 
(19) B 2 0 = 15.0 


Table 2: Frequency of Oscillation (KHz) 



M=4.18 

M=4.48 

M=4.79 

Experiment 13 

148 

425 

720 

Present method 

163 

431 

701-716 

Matsuo et a/. 4 

160 

- 

725 

Wilson & Sussman 3 

- 

- 

530 

Hosangadi et a/. 23 

- 

- 

450-500 
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Reaction Shocks 



(a) 



Rarefaction Wave 


(b) 


Reaction Front 
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(C) 


Figure 1 . Schematic diagram of one dimensional 
wave interactions: (a) generation of reaction shocks 
and interaction with bow shock; (b) initiation of new 
reaction front, and extinguishment of original reaction 
front: (c) interaction of rarefaction wave with the bow 
shock. 


Figure 2. McVey and Tong 16 x-t diagram along the 
stagnation streamline. 
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Density 
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(b) 


Pressure 


ICO Mass Fraction 




(c) 


(d) 
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Figure 3. Experimental and computational results for a projectile moving at M = 4.18 in a stoichiometric 
mixture of hydrogen-air: (a) experimental shadowgraph image (Lehr 15 ); (b) density (p//>oo), contour range: 
(min=.25, max=5.35, inc=.l); (c) pressure (p/poo), (1.5, 23.0, .25); (d) H 2 0 mass fraction, (.005, .22, .005). 
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Density 


Pressure 




t = 13.57 iisec 


t = 0 


(a) 


(b) 


Figure 4. History of density and pressure along the stagnation streamline; M = 4.18: (a) density (p/poo), 
contour range: (min=2.5, max=5.85, inc=.05); (b) pressure (p/poo), (18.0, 28.1, .25). 



(a) 


Density 



Figure 5. Experimental and computational results for a projectile moving at M = 4.48 in a stoichiometric 
mixture of hydrogen-air: (a) experimental shadowgraph image (Lehr 15 ) taken at a 12° angle from the 
perpendicular direction; (b) density (p/poo), contour range: (min=.35, max=6.35, inc=.l). 
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Pressure 



H20 Mass Fraction 



Figure 7. Experimental and computational results for a projectile moving at M = 4.79 in a stoichiometric 
mixture of hydrogen-air: (a) experimental shadowgraph image (Lehr 15 ); (b) density (p/poo), contour range: 
(min=.35, max=5.65, inc=.l); (c) pressure (p/poo), (-75, 30.75, .5); (d) H 2 O mass fraction, (.005, .225, .005). 
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